# This script reproduces Table 6 #
rm(list=ls())
options(digits=4)
pkg <- c("PanelMatch", "ggplot2", "plm", 
         "matrixStats",
         "foreign", 
         "stringi", "stringdist",
         "DataCombine","lmtest", "multiwayvcov",
         "data.table",
         "dplyr", 
         "stargazer",
         "clubSandwich", "mediation",
         "doBy")
lapply(pkg, require, character.only = TRUE)

# load datasets
proprietary <- read.csv("proprietary.csv")
load("non_p.RData")

# merge 
d_sub1 <- merge(d_sub1, proprietary, by = c("county_id", "time"), all.x = T)

# with lag
model_area_con <- lm(log_success_area_f1 ~ 
                       tour + 
                       log_success_area_l1 +
                       log_success_area_l2 +
                       log_success_area_l3 +
                       log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                       log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                       log_expd_l1 + log_middle_school_l1 + 
                       msec2currentsec_l1 + 
                       msec2currentgvn_l1 + 
                       mayor2currentsec_l1 + 
                       mayor2currentgvn_l1 + as.factor(county_id) + as.factor(time),
                     data = d_sub1)
vcov_area_con <- clubSandwich::vcovCR(model_area_con, type="CR1S", cluster = d_sub1$county_id)
se_area_con <- coeftest(model_area_con, vcov_area_con)
reg_results_lag_land <- list("model" = model_area_con, 
                             "se" = se_area_con)

# without lag
model_area_con <- lm(log_success_area_f1 ~ 
                       tour + 
                       log_gdppc_l1 + log_gdp_l1 + log_fia_l1 + log_loan_l1 +
                       log_hos_l1 + log_industrial_l1 + log_rev_l1 +
                       log_expd_l1 + log_middle_school_l1 + 
                       msec2currentsec_l1 + 
                       msec2currentgvn_l1 + 
                       mayor2currentsec_l1 + 
                       mayor2currentgvn_l1 + as.factor(county_id) + as.factor(time),
                     data = d_sub1)
vcov_area_con <- clubSandwich::vcovCR(model_area_con, type="CR1S", cluster = d_sub1$county_id)
se_area_con <- coeftest(model_area_con, vcov_area_con)
reg_results_land <- list("model" = model_area_con, 
                         "se" = se_area_con)


sink("output/Table6.tex")
stargazer(reg_results_land$model,
          reg_results_lag_land$model,
          font.size = "large", digits = 4, 
          type = "latex",
          dep.var.labels = c("Logged Successfully Transacted Area Next Month",
                             "Logged Successfully Transacted Area Next Month"),
          column.labels = c("Logged Successfully Transacted Area Next Month",
                            "Logged Successfully Transacted Area Next Month"),
          dep.var.labels.include = F,
          style = "qje", no.space = FALSE, keep.stat = c("n", "rsq"),
          omit = c("county_id", "time"
          ),
          omit.yes.no = c("No", "Yes"),
          add.lines = list(c("County FE?", "Yes", "Yes"),
                           c("Time FE?", "Yes", "Yes"),
                           c("Controls?", "Yes", "Yes"),
                           c("Lagged dependent variable included?", "No", "Yes")),
          covariate.labels = c("Inspection",
                               "Logged successfully transacted area last month",
                               "Logged successfully transacted area two months prior",
                               "Logged successfully transacted area three months prior",
                               "logged GDP per capita last year", 
                               "logged total GDP last year", 
                               "logged fixed asset investment last year",
                               "logged volume of loans last year", 
                               "logged number of hospitals last year", 
                               "logged amount of industrial output last year", 
                               "logged revenue last year", 
                               "logged expenditure last year", 
                               "logged number of secondary schools last year",
                               "political connections 1 last year",
                               "political connections 2 last year",
                               "political connections 3 last year",
                               "political connections 4 last year"),
          star.cutoffs = c(.10, .05, .01), 
          se = list(reg_results_land$se[,2], 
                    reg_results_lag_land$se[,2]),
          p = list(reg_results_land$se[,4], 
                   reg_results_lag_land$se[,4]),
          header = F, notes = "robust standard errors clustered by county in parentheses", 
          float = F)
sink()

